function pr=prob_ReSTAT_VARvcsav(gg_cub,bd_cub,bdtype_cub,fchs_cub,ind_cub,idx_sw_cub,curdom_cub,mktcat_cub,hhi_cub,pop65_cub,eerr_cub,df_cub,pre,act,theta)   %For 26 parameters with loop  %###UPDATED on 02/09/2016 switching cost expressed as positive
%### UPDATED on 03/17/2017: include the ln(pr) as the expected error for
%future flows.
%### UPDATED on 05/05/2017: to measure mktdom and swc, use the quantile of
%bed size, which are mktdom*bd[1-4] and swc*bd[1-4]; include "bdtype" in
%the inputs.
%### UPDATED on 07/16/2018: replace mktdom dummy with market share.
%### UPDATED on 04/23/2020: this is the same specification but using new
%simulated future states that include obs and unobs market characteristics.
%### UPDATED on 05/05/2020: introduce decreasing sunk cost VARYing by
%vendor
%### UPDATED on 06/22/2020: simplify by reducing for loop
%### UPDATED on 11/01/2020: use concurrent mktdom for variable cost savings

ev1=gg_cub.*theta(end-bdtype_cub-3);
ev1(bd_cub==0)=0;
% cost-saving per bed
vbd=theta(fchs_cub);
vbd(fchs_cub==11)=0;
vbd(fchs_cub==12)=0;
ev2=bd_cub.*vbd;

% variable cost savings X mktdom %### UPDATED on 12/19/2020: VARY by
% vendors; Vendor 2 and 13 have NO coefficients here
%ev_vcdom=bd_cub.*curdom_cub*theta(52);
coef_VARvcsav=theta(fchs_cub+38);
coef_VARvcsav(fchs_cub==11|fchs_cub==12|fchs_cub==13|fchs_cub==24)=0;
ev_vcdom=bd_cub.*curdom_cub.*coef_VARvcsav;

% sunk cost
ev3=zeros(size(ind_cub));
ev3(ind_cub>0)=theta(ind_cub(ind_cub>0)-1);
%vdecsunk=theta(ind_cub+23);   %### UPDATED on 05/05/2020: updated sunk cost to include the decreasing part VARYing by vendor 
%vdecsunk(ind_cub==0)=0;
%vdecsunk(ind_cub==1)=0;
%ev3=ev3_part1+vdecsunk.*(sunk2>0).*sunk2+MKTDOM_cub.*sunk2*theta(37);
% switching cost
ev4=idx_sw_cub.*theta(end+bdtype_cub-4);
ev4(bd_cub==0)=0;
%market FEs; coefficients for market categories: theta(49, 50, 51)
%ev_mktcat=theta(mktcat_cub+24);
ev_mktcat=theta(mktcat_cub+48);  % UPDATED on 07/04/2020: the position of the parameters changes
ev_mktcat(fchs_cub==11|fchs_cub==12|mktcat_cub==0)=0;  
% HHI, varying by vendors
vhhi=theta(fchs_cub+12);
vhhi(fchs_cub==11|fchs_cub==12)=0;
ev_hhi=hhi_cub.*vhhi;
% totpop65, varying by vendors
vpop65=theta(fchs_cub+24);
vpop65(fchs_cub==11|fchs_cub==12)=0;
ev_pop65=pop65_cub.*vpop65;
% SUM up
%pv=sum((ev1+ev2+ev3+ev4+eerr_cub).*df_cub,2); 
eev=mean(sum((ev1+ev2+ev3+ev4+ev_mktcat+ev_hhi+ev_pop65+ev_vcdom-eerr_cub).*df_cub,2),3);
v=reshape(eev,13,[])';
expv=exp(v);
expv(pre(:,1)>1,1)=0;
idx=sub2ind(size(expv),(1:size(expv,1))',act);
pr=expv(idx)./sum(expv,2);
ll=-sum(log(pr));





